Tutorial Part IV - Movement - Solution Exercise 2
Exercise 4
Exercise 4.3
Three animals have been collared and tracked, and we want to know their home range size and whether the animals actually avoid agriculture and prefer forest.
There is telemetry data available from three animals (DieTiere_32633.shp) in the folder data-raw. Additionally, land cover (Landnutzung_3035.shp) and rivers (linearegew_32633.shp) are available in the geo-raw folder. Land use has a column NS1 with coded info on 1 (urban areas), 2 (agriculture), 3 (forest), 4 (swamps), and 5 (water bodies).
After loading the necessary libraries, the data and making a plot, please do the following:
- Calculate their home ranges.
- Calculate their home ranges.
- Buffer the home range with the mean displacement per step of all three animals.
- Buffer the home range with the mean displacement per step of all three animals.
- Clip the home ranges with the underlying landscape and calculate the percentage of land use classes 2 (agriculture) and 2 (forest) within the home range (hint: use
st_intersection()).
- Clip the home ranges with the underlying landscape and calculate the percentage of land use classes 2 (agriculture) and 2 (forest) within the home range (hint: use
- In which land use are the animals actually located? Extract land use underneath each telemetry location (hint:
st_join()). Do the animals actually avoid agriculture?
- In which land use are the animals actually located? Extract land use underneath each telemetry location (hint:
Get started
## Open libraries
## general basics
library(here)
library(ggplot2)
## the basic packages for spatial analyses and viz
library(sf)
library(terra)
library(tmap)
## for telemetry and home range analyses
library(move)
library(adehabitatHR)
library(sp)Define workspace and folder structure:
work_wd <- here()
data_raw_wd <- here("data-raw")
geo_raw_wd <- here("data-raw", "geo-raw")
output_wd <- here("output")Load data sets:
# if you are not sure about the file names, check the folder content with
# list.files(data_raw_wd)
animals <- st_read(layer = "DieTiere_32633", dsn = data_raw_wd) ## Reading layer `DieTiere_32633' from data source
## `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\data-raw'
## using driver `ESRI Shapefile'
## Simple feature collection with 388 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 390397.5 ymin: 5862042 xmax: 396748.2 ymax: 5868366
## Projected CRS: WGS 84 / UTM zone 33N
landcov <- st_read(layer = "Landnutzung_3035", dsn = geo_raw_wd) ## Reading layer `Landnutzung_3035' from data source
## `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\data-raw\geo-raw'
## using driver `ESRI Shapefile'
## Simple feature collection with 182 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4539362 ymin: 3304861 xmax: 4566087 ymax: 3324860
## Projected CRS: ETRS89-extended / LAEA Europe
rivers <- st_read(layer = "linearegew_32633", dsn = geo_raw_wd) ## Reading layer `linearegew_32633' from data source
## `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\data-raw\geo-raw'
## using driver `ESRI Shapefile'
## Simple feature collection with 2834 features and 13 fields
## Geometry type: MULTILINESTRING
## Dimension: XYZ
## Bounding box: xmin: 381989.2 ymin: 5852103 xmax: 408856.6 ymax: 5871804
## z_range: zmin: 0 zmax: 0
## Projected CRS: WGS 84 / UTM zone 33N
Define or transform coordinate reference system (CRS):
lc_32633 <- st_transform(landcov, 32633)Make a plot of the data:
## Create a color code for the 5 land use classes
land_cols <- c("red", "antiquewhite", "chartreuse4", "cornsilk4", "blue")
## TODO - argh, Punkte und Linien erscheinen nicht im html....
## CED: Hab leider keinen Schimmer von base pltos... ggplot for the rescue?
plot(lc_32633[,"NS1"], col = land_cols[lc_32633$NS1], border = "transparent")
plot(animals[,1], col = animals$ID, add = TRUE)
plot(rivers[,1], col = 'blue', lwd = 0.1, add = TRUE)ggplot2:
## Create a color code for the 5 land use classes
land_cols <- c("red", "antiquewhite", "chartreuse4", "cornsilk4", "navyblue")
## Create a color code for the 3 individuals
id_cols <- viridis::inferno(n = 12)[c(2,4,6)]
ggplot() +
geom_sf(data = lc_32633, aes(fill = as.factor(NS1)), color = NA, alpha = .4) +
geom_sf(data = rivers, color = "navyblue", alpha = .4) +
geom_sf(data = animals, aes(color = as.factor(ID))) +
scale_color_manual(values = id_cols, name = "Animal ID") +
scale_fill_manual(values = land_cols, name = "Land Use Class") +
theme_minimal(base_size = 15)tmap:
## Create a color code for the 5 land use classes
land_cols <- c("red", "antiquewhite", "chartreuse4", "cornsilk4", "navyblue")
## Create a color code for the 3 individuals
id_cols <- viridis::inferno(n = 12)[c(2,4,6)]
tmap_mode(mode = "view")
tm_shape(shp = lc_32633[,'NS1']) +
tm_polygons(col = "NS1", palette = land_cols, alpha = .4) +
tm_shape(shp = rivers[,1]) +
tm_lines(col = 'blue') +
tm_shape(shp = animals[,'ID']) +
tm_dots(col = "ID", palette = id_cols)- a. Calculate their home ranges.
and - b. Buffer the home range with the mean displacement per step of all three animals.
First , we retrieve the step length. For this, transform the sf object into a move object:
anim_coords <- as.data.frame(st_coordinates(animals)) ## returns a list
animals$X <- anim_coords$X
animals$Y <- anim_coords$Y
animals_ng <- st_drop_geometry(animals)
animals_move <- move(x = animals_ng$X, y = animals_ng$Y,
time = as.POSIXct(animals_ng$MYDATUM,
format = "%Y-%m-%d 8:0:0", tz = "CET"), ## CED: Puh, das ist aber komplex... geht's nicht iwie einfacher, z.B. mit lubridate?
proj = CRS("+init=epsg:32633"),
data = animals_ng, animal = animals_ng$ID, sensor = "GPS")Check the data by plotting the new move object:
plot(animals_move, type = "l")Since we have 3 animals, a list with 3 elements is returned:
steplength_per_animal <- distance(animals_move)
mean_stepl_p_anim <- lapply(steplength_per_animal,mean)
(mean(unlist(mean_stepl_p_anim)))## [1] 303.0478
The mean step length is approx. 300 m.
Now use this value to buffer the MCPs:
## for using mcp() in {adehabitatHR}, convert to SpatialPointsDataframe of {sp}
animals_sp <- as(animals, "Spatial")
## calculate MCP
mcp <- mcp(animals_sp[,"ID"], percent = 95)
## the buffer command needs however an sf object
mcp_sf <- as(mcp, "sf")
mcp_buf <- st_buffer(mcp_sf, 300)
## the plot is correct
plot(mcp)
plot(animals_sp, add = TRUE, col = animals_sp$ID)
plot(mcp_buf[,1], add = TRUE)
plot(animals_sp, add = TRUE, col = id_cols)
plot(mcp, add = TRUE)## save buffer to disc
st_write(obj = mcp_buf,
dsn = output_wd,
layer = "mcp_buffer_32633",
driver = "ESRI Shapefile",
delete_layer = TRUE)## Deleting layer `mcp_buffer_32633' using driver `ESRI Shapefile'
## Writing layer `mcp_buffer_32633' to data source
## `C:/Users/DataVizard/Google Drive/Work/Eco/Projects/2022_BioDivCourse/Course4_MoveQ/output' using driver `ESRI Shapefile'
## Writing 3 features with 2 fields and geometry type Polygon.
- c. Clip the home ranges with the underlying landscape and calculate the percentage of land use classes 2 (agriculture) and 2 (forest) within the home range (hint: use
st_intersection()).
## check the total area of the buffer in km2 per animal:
buf_area <- st_area(mcp_buf) / 1e6
## clip with land use
mcp_buf_lc <- st_intersection(lc_32633, mcp_buf)
## have a look at the object:
mcp_buf_lc## Simple feature collection with 45 features and 14 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 390158.2 ymin: 5861742 xmax: 397048.2 ymax: 5868666
## Projected CRS: WGS 84 / UTM zone 33N
## First 10 features:
## AREA PERIMETER CLCGES_BES CLCGES_B_1 NS BRD SZENE SPHEROID TKNR NS1
## 2 33397200 79888.10 23163 23162 312 1 19323 Krassovsky N33111 3
## 5 11562800 44001.70 23265 23264 211 1 19323 Krassovsky N33111 2
## 15 18208500 73976.60 23764 23763 231 1 19323 Krassovsky N33111 2
## 16 8568800 37160.80 23788 23787 231 1 19323 Krassovsky N33111 2
## 28 5142930 13701.10 24028 24027 312 1 19323 Krassovsky N33111 3
## 35 1227710 4956.00 24314 24313 243 1 19323 Krassovsky N33111 2
## 36 1448920 11624.70 24348 24347 311 1 19323 Krassovsky N33111 3
## 41 654073 3560.28 24431 24430 211 1 19323 Krassovsky N33111 2
## 42 3534710 11425.00 24454 24453 311 1 19323 Krassovsky N33111 3
## 44 46474200 125933.00 24480 24479 211 1 19323 Krassovsky N33111 2
## NS2 NS3 id area geometry
## 2 31 312 1 1444.323 MULTIPOLYGON (((394196.9 58...
## 5 21 211 1 1444.323 POLYGON ((392577.5 5866547,...
## 15 23 231 1 1444.323 POLYGON ((391413.6 5866716,...
## 16 23 231 1 1444.323 POLYGON ((395139.3 5866196,...
## 28 31 312 1 1444.323 POLYGON ((390833 5866283, 3...
## 35 24 243 1 1444.323 POLYGON ((392292.4 5866469,...
## 36 31 311 1 1444.323 POLYGON ((395809.1 5865790,...
## 41 21 211 1 1444.323 POLYGON ((394729.2 5865833,...
## 42 31 311 1 1444.323 POLYGON ((390395.2 5865319,...
## 44 21 211 1 1444.323 MULTIPOLYGON (((392247.2 58...
## take care, you have to recalculate the area!
mcp_buf_lc$areakm2 <- st_area(mcp_buf_lc) / 1e6
## now, for each animal id, summarize the area per NS1-category 2 and 3:
agri <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 2], ## CED: Das können wir wohl von denen nicht erwarten oder??
mcp_buf_lc$id[mcp_buf_lc$NS1 == 2],
FUN = sum)
forest <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 3],
mcp_buf_lc$id[mcp_buf_lc$NS1 == 3],
FUN = sum)
## percentages
round(agri / buf_area * 100, 0)## Units: [1/m^2]
## 1 2 3
## 71 48 85
round(forest / buf_area * 100, 0)## Units: [1/m^2]
## 1 2 3
## 29 52 14
- d. In which land use are the animals actually located? Extract land use underneath each telemetry location (hint:
st_join()). Do the animals actually avoid agriculture?
spatjoin <- st_join(animals, lc_32633)
myres <- table(spatjoin$NS1, spatjoin$ID)
myres##
## 1 2 3
## 2 48 25 95
## 3 147 64 9
## now divide that by the total number of locations per animal
myn <- table(animals$ID)
round((myres[1,] / myn*100), 0) ## agriculture##
## 1 2 3
## 25 28 91
round((myres[2,] / myn*100), 0) ## forest##
## 1 2 3
## 75 72 9
Compare the values with the percentages above.
For example, animal 1 had 71% of agriculture and 29% of forest in its HR, however, only 25% of the locations were actually in agriculture, but 75% of its locations were in forest. We could conclude that animal 1 avoids agriculture.
Session Info
Sys.time()## [1] "2022-03-07 18:21:20 CET"
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
## [5] LC_TIME=C
## system code page: 65001
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] adehabitatHR_0.4.19 adehabitatLT_0.3.25 CircStats_0.2-6
## [4] boot_1.3-28 MASS_7.3-54 adehabitatMA_0.3.14
## [7] ade4_1.7-18 deldir_1.0-6 move_4.1.6
## [10] rgdal_1.5-28 raster_3.5-11 sp_1.4-6
## [13] geosphere_1.5-14 tmap_3.3-2 terra_1.4-22
## [16] sf_1.0-5 ggplot2_3.3.5 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-2 httr_1.4.2 rprojroot_2.0.2
## [4] tools_4.1.2 bslib_0.3.1 utf8_1.2.2
## [7] R6_2.5.1 KernSmooth_2.23-20 DBI_1.1.2
## [10] colorspace_2.0-2 withr_2.4.3 gridExtra_2.3
## [13] tidyselect_1.1.1 leaflet_2.0.4.1 compiler_4.1.2
## [16] leafem_0.1.6 textshaping_0.3.6 xml2_1.3.2
## [19] bookdown_0.24 sass_0.4.0 scales_1.1.1
## [22] classInt_0.4-3 proxy_0.4-26 systemfonts_1.0.3
## [25] stringr_1.4.0 digest_0.6.29 rmarkdown_2.11
## [28] base64enc_0.1-3 dichromat_2.0-0 pkgconfig_2.0.3
## [31] htmltools_0.5.2 highr_0.9 fastmap_1.1.0
## [34] htmlwidgets_1.5.4 rlang_0.4.12 farver_2.1.0
## [37] jquerylib_0.1.4 generics_0.1.1 jsonlite_1.7.2
## [40] crosstalk_1.2.0 dplyr_1.0.7 magrittr_2.0.1
## [43] s2_1.0.7 Rcpp_1.0.7 munsell_0.5.0
## [46] fansi_0.5.0 viridis_0.6.2 abind_1.4-5
## [49] lifecycle_1.0.1 stringi_1.7.5 leafsync_0.1.0
## [52] yaml_2.2.1 tmaptools_3.1-1 grid_4.1.2
## [55] parallel_4.1.2 crayon_1.4.2 lattice_0.20-45
## [58] stars_0.5-5 knitr_1.36 pillar_1.6.4
## [61] codetools_0.2-18 wk_0.6.0 XML_3.99-0.8
## [64] glue_1.4.2 evaluate_0.14 leaflet.providers_1.9.0
## [67] png_0.1-7 vctrs_0.3.8 rmdformats_1.0.3
## [70] gtable_0.3.0 purrr_0.3.4 assertthat_0.2.1
## [73] cachem_1.0.6 xfun_0.27 lwgeom_0.2-8
## [76] e1071_1.7-9 ragg_1.1.3 class_7.3-19
## [79] viridisLite_0.4.0 tibble_3.1.6 memoise_2.0.1
## [82] units_0.7-2 ellipsis_0.3.2